Licensed under the Creative Commons attribution-noncommercial license, http://creativecommons.org/licenses/by-nc/3.0/. Please share and remix noncommercially, mentioning its origin.
CC-BY_NC Produced in R version 3.2.1 using pomp version 0.72.2.

In this example, we revisit the analysis of He et al. (2010). This is intended to be a concrete example of a full data analysis. Whilst He et al. (2010) examined measles in 20 cities, we will focus on London alone.

Measles revisited

Challenges & motivation

  • Understanding, forecasting, managing epidemiological systems increasingly depends on models
  • Increasing interest in opening up the black box: mechanistic models
  • Real epidemiological systems:
    • are nonlinear
    • are stochastic
    • are nonstationary
    • evolve in continuous time
    • have hidden variables
    • can be measured only with (large) error
  • Dynamics of infectious disease outbreaks illustrate this well

Outline

  • revisit classic measles data set
  • use a model that
    1. expresses our current understanding of measles dynamics
    2. cannot be fit by existing likelihood-based methods
  • examine data from large and small towns using the same model
  • does our perspective on this disease change?
  • what bigger lessons can we learn regarding inference for dynamical systems?

He et al. (2010)

Data sets

  • Twenty towns
    • population sizes: 2k–3.4M
    • Weekly case reports, 1950–1963
    • Annual birth records and population sizes, 1944–1963

Continuous-time Markov process model

  • \(B(t) = \text{birth rate, from data}\)
  • \(N(t) = \text{population size, from data}\)

  • \(\mathrm{cases}_t\,\vert\,I{\to}R=z_t \sim \mathrm{normal}\left(\rho\,z_t,\rho\,(1-\rho)\,z_t+(\psi\,\rho\,z_t)^2\right)\)

  • entry into susceptible class: \[\mu_{BS}(t) = (1-c)\,B(t-\tau)+c\,\delta(t-t_0)\,\int_{t-1}^{t}\,B(t-\tau-s)\,ds\]

  • force of infection \[\mu_{SE}(t) = \tfrac{\beta(t)}{N(t)}\,(I+\iota)\,\zeta(t)\]

  • \(c = \text{cohort effect}\)
  • \(\tau = \text{school-entry delay}\)
  • \(\beta(t) = \text{school-term transmission}\)
  • \(\iota = \text{imported infections}\)

  • \(\zeta(t) = \text{white noise with intensity}\,\sigma_{SE}\)

Fitting procedure

  • a large Latin-hypercube design was used to initiate searches
  • iterated filtering to maximize the likelihood
  • point estimates of all parameters for 20 cities
  • profile likelihoods to quantify uncertainty in London and Hastings

Imported infections

\[\mu_{SE}=\frac{\beta(t)}{N(t)}\,(I+\iota)\,\zeta(t)\]

Seasonality

Notable findings

Parameter estimates

pop R0 amplitude LP IP alpha iota rho psi sigmaSE
Halesworth 2170 33.100 0.381 7.870 2.280 0.948 0.00912 0.754 0.641 0.0748
Lees 4250 29.700 0.153 8.510 2.050 0.968 0.03110 0.612 0.681 0.0802
Mold 6410 21.400 0.271 5.930 1.780 1.040 0.01450 0.131 2.870 0.0544
Dalton.in.Furness 10600 28.300 0.203 5.480 1.980 0.989 0.03860 0.455 0.818 0.0779
Oswestry 11000 52.900 0.339 10.300 2.720 1.040 0.02980 0.631 0.476 0.0699
Northwich 18300 30.100 0.423 8.510 3.010 0.948 0.06020 0.795 0.402 0.0857
Bedwellty 28900 24.700 0.160 6.820 3.030 0.937 0.03960 0.311 0.951 0.0611
Consett 39100 35.900 0.200 9.070 2.660 1.010 0.07310 0.650 0.406 0.0712
Hastings 65700 34.200 0.299 7.000 5.440 1.000 0.18600 0.695 0.396 0.0955
Cardiff 245000 34.400 0.223 9.860 3.090 0.996 0.14100 0.602 0.270 0.0539
Bradford 294000 32.100 0.236 8.510 3.360 0.991 0.24400 0.599 0.190 0.0451
Hull 302000 38.900 0.221 9.180 5.460 0.968 0.14200 0.582 0.256 0.0636
Nottingham 307000 22.600 0.157 5.720 3.690 0.982 0.17000 0.609 0.258 0.0380
Bristol 443000 26.800 0.203 6.190 4.940 1.010 0.44100 0.626 0.201 0.0392
Leeds 510000 47.800 0.267 9.480 10.900 1.000 1.25000 0.666 0.167 0.0778
Sheffield 515000 33.100 0.313 7.230 6.380 1.020 0.85300 0.649 0.175 0.0428
Manchester 704000 32.900 0.290 11.100 6.940 0.965 0.59000 0.550 0.161 0.0551
Liverpool 802000 48.100 0.305 7.900 9.800 0.978 0.26300 0.494 0.136 0.0533
Birmingham 1120000 43.400 0.428 8.520 11.600 1.010 0.34300 0.544 0.178 0.0611
London 3390000 56.800 0.554 13.100 12.500 0.976 2.90000 0.488 0.116 0.0878
r NA 0.444 0.293 0.349 0.946 0.108 0.93200 -0.203 -0.931 -0.3170

\(r = \mathrm{cor}(\log{\hat\theta},\log{N_{1950}})\)

Extrademographic stochasticity

\[\mu_{SE}=\frac{\beta(t)}{N(t)}\,(I+\iota)\,\zeta(t)\]

Reporting rate

Cohort effect

\(R_0\)

Birth delay

Predicted vs observed critical community size

Practicum

Preliminaries

To get started, we must install pomp if it is not already installed. The following commands install pomp from source (on github).

require(devtools)
install_github("kingaa/pomp")

If we were to run into trouble in the above, we would consult the package website.

Now we’ll load the packages we’ll need, and set the random seed, to allow reproducibility.

set.seed(594709947L)
require(ggplot2)
theme_set(theme_bw())
require(plyr)
require(reshape2)
require(magrittr)
require(pomp)
stopifnot(packageVersion("pomp")>="0.70-1")

Data and covariates

Now we’ll load the data and covariates. The data are measles reports from 20 cities in England and Wales. We also have information on the population sizes and birth-rates in these cities; we’ll treat these variables as covariates.

baseurl <- "http://kinglab.eeb.lsa.umich.edu/SBIED/"
download.file(paste0(baseurl,"data/twentycities.rda"),destfile="./twentycities.rda")
load("twentycities.rda")
measles %>% 
  mutate(year=as.integer(format(date,"%Y"))) %>%
  subset(town=="London" & year>=1950 & year<1964) %>%
  mutate(time=(julian(date,origin=as.Date("1950-01-01")))/365.25+1950) %>%
  subset(time>1950 & time<1964, select=c(time,cases)) -> dat
demog %>% subset(town=="London",select=-town) -> demog

Let’s plot the data and covariates.

dat %>% ggplot(aes(x=time,y=cases))+geom_line()

demog %>% melt(id="year") %>%
  ggplot(aes(x=year,y=value))+geom_point()+
  facet_wrap(~variable,ncol=1,scales="free_y")

demog %>% 
  summarize(
    time=seq(from=min(year),to=max(year),by=1/12),
    pop=predict(smooth.spline(x=year,y=pop),x=time)$y,
    birthrate=predict(smooth.spline(x=year+0.5,y=births),x=time-4)$y
    ) -> covar
plot(pop~time,data=covar,type='l')
points(pop~year,data=demog)

plot(birthrate~time,data=covar,type='l')
points(births~year,data=demog)

plot(birthrate~I(time-4),data=covar,type='l')
points(births~I(year+0.5),data=demog)

The partially observed Markov process model

The (unobserved) process model

Let’s evaluate the hypothesis that these data were generated by an SEIR model. The SEIR model is a compartmental model that, diagrammatically, looks as follows.

\(b = \text{births}\)
\(S = \text{susceptibles}\)
\(E = \text{exposed, incubating}\)
\(I = \text{infectious}\)
\(R = \text{recovered}\)

We require a simulator for this model. The following implements a simulator.

rproc <- Csnippet("
  double beta, br, seas, foi, dw, births;
  double rate[6], trans[6];
  
  // cohort effect
  if (fabs(t-floor(t)-251.0/365.0) < 0.5*dt) 
    br = cohort*birthrate/dt + (1-cohort)*birthrate;
  else 
    br = (1.0-cohort)*birthrate;

  // term-time seasonality
  t = (t-floor(t))*365.25;
  if ((t>=7&&t<=100) || (t>=115&&t<=199) || (t>=252&&t<=300) || (t>=308&&t<=356))
      seas = 1.0+amplitude*0.2411/0.7589;
    else
      seas = 1.0-amplitude;

  // transmission rate
  beta = R0*seas*(1.0-exp(-(gamma+mu)*dt))/dt;
  // force of infection
  foi = beta*pow(I+iota,alpha)/pop;
  // white noise (extrademographic stochasticity)
  dw = rgammawn(sigmaSE,dt);

  rate[0] = dw*foi/dt;  //infection rate (stochastic)
  rate[1] = mu;             // natural S death
  rate[2] = sigma;        // rate of ending of latent stage
  rate[3] = mu;             // natural E death
  rate[4] = gamma;        // recovery
  rate[5] = mu;             // natural I death

  // Poisson births
  births = rpois(br*dt);
  
  // transitions between classes
  reulermultinom(2,S,&rate[0],dt,&trans[0]);
  reulermultinom(2,E,&rate[2],dt,&trans[2]);
  reulermultinom(2,I,&rate[4],dt,&trans[4]);

  S += births   - trans[0] - trans[1];
  E += trans[0] - trans[2] - trans[3];
  I += trans[2] - trans[4] - trans[5];
  R = pop - S - E - I;
  W += (dw - dt)/sigmaSE;  // standardized i.i.d. white noise
  C += trans[4];           // true incidence
")

In the above, \(C\) represents the true incidence, i.e., the number of new infections occurring over an interval. Since recognized measles infections are quarantined, we argue that most infection occurs before case recognition so that true incidence is a measure of the number of individuals progressing from the I to the R compartment in a given interval.

We complete the process model definition by specifying the distribution of initial unobserved states. The following codes assume that the fraction of the population in each of the four compartments is known.

initlz <- Csnippet("
  double m = pop/(S_0+E_0+I_0+R_0);
  S = nearbyint(m*S_0);
  E = nearbyint(m*E_0);
  I = nearbyint(m*I_0);
  R = nearbyint(m*R_0);
  W = 0;
  C = 0;
")

The measurement model

We’ll model both under-reporting and measurement error. We want \(\mathbb{E}[\text{cases}|C] = \rho\,C\), where \(C\) is the true incidence and \(0<\rho<1\) is the reporting efficiency. We’ll also assume that \(\mathrm{Var}[\text{cases}|C] = \rho\,(1-\rho)\,C + (\psi\,\rho\,C)^2\), where \(\psi\) quantifies overdispersion. Note that when \(\psi=0\), the variance-mean relation is that of the binomial distribution. To be specific, we’ll choose \(\text{cases|C} \sim f(\cdot|\rho,\psi,C)\), where \[f(c|\rho,\psi,C) = \Phi(c+\tfrac{1}{2},\rho\,C,\rho\,(1-\rho)\,C+(\psi\,\rho\,C)^2)-\Phi(c-\tfrac{1}{2},\rho\,C,\rho\,(1-\rho)\,C+(\psi\,\rho\,C)^2),\] where \(\Phi(x,\mu,\sigma^2)\) is the c.d.f. of the normal distribution with mean \(\mu\) and variance \(\sigma^2\).

The following computes \(\mathbb{P}[\text{cases}|C]\).

dmeas <- Csnippet("
  double m = rho*C;
  double v = m*(1.0-rho+psi*psi*m);
  double tol = 1.0e-18;
  if (cases > 0.0) {
      lik = pnorm(cases+0.5,m,sqrt(v)+tol,1,0)-pnorm(cases-0.5,m,sqrt(v)+tol,1,0)+tol;
  } else {
    lik = pnorm(cases+0.5,m,sqrt(v)+tol,1,0)+tol;
  }
")

The following codes simulate \(\text{cases} | C\).

rmeas <- Csnippet("
  double m = rho*C;
  double v = m*(1.0-rho+psi*psi*m);
  double tol = 1.0e-18;
  cases = rnorm(m,sqrt(v)+tol);
  if (cases > 0.0) {
    cases = nearbyint(cases);
  } else {
    cases = 0.0;
  }
")

Constructing the pomp object

We put all the model components together with the data in a call to pomp:

dat %>% 
  pomp(t0=with(dat,2*time[1]-time[2]),
       time="time",
       rprocess=euler.sim(rproc,delta.t=1/365.25),
       initializer=initlz,
       dmeasure=dmeas,
       rmeasure=rmeas,
       covar=covar,
       tcovar="time",
       zeronames=c("C","W"),
       statenames=c("S","E","I","R","C","W"),
       paramnames=c("R0","mu","sigma","gamma","alpha","iota",
                    "rho","sigmaSE","psi","cohort","amplitude",
                    "S_0","E_0","I_0","R_0")
       ) -> m1

The following codes plot the data and covariates together.

m1 %>% as.data.frame() %>% 
  melt(id="time") %>%
  ggplot(aes(x=time,y=value))+
  geom_line()+
  facet_grid(variable~.,scales="free_y")

He et al. (2010) estimated the parameters of this model. The following attaches the

readRDS("He2010_mles.rds") %>% 
  subset(town=="London") -> mle
paramnames <- c("R0","mu","sigma","gamma","alpha","iota",
                "rho","sigmaSE","psi","cohort","amplitude",
                "S_0","E_0","I_0","R_0")
theta <- unlist(mle[paramnames])
kable(t(mle))
23
town London
loglik -3804.935
loglik.sd 0.1570628
mu 0.02
delay 4
sigma 28.88585
gamma 30.40555
rho 0.487543
R0 56.78228
amplitude 0.5540616
alpha 0.9755556
iota 2.895151
cohort 0.5570074
psi 0.1158933
S_0 0.02970213
E_0 5.17412e-05
I_0 5.1426e-05
R_0 0.9701947
sigmaSE 0.0877794

Verify that we get the same likelihood as He et al. (2010).

require(foreach)
require(doMC)
registerDoMC()

set.seed(998468235L,kind="L'Ecuyer")
mcopts <- list(preschedule=FALSE,set.seed=TRUE)
paropts <- list(.options.multicore=mcopts)

foreach(i=1:4,
        .packages="pomp",
        .options.multicore=mcopts) %dopar% {
  pfilter(m1,Np=10000,params=theta)
} -> pfs
logmeanexp(sapply(pfs,logLik),se=TRUE)
##                          se 
## -3801.8069447     0.7678851

Simulations at the MLE.

m1 %>% 
  simulate(params=theta,nsim=10,as.data.frame=TRUE,include.data=TRUE) %>%
  ggplot(aes(x=time,y=cases,group=sim,color=(sim=="data")))+
  guides(color=FALSE)+
  geom_line()+facet_wrap(~sim,ncol=2)

m1 %>% 
  simulate(params=theta,nsim=100,as.data.frame=TRUE,include.data=TRUE) %>%
  subset(select=c(time,sim,cases)) %>%
  mutate(data=sim=="data") %>%
  ddply(~time+data,summarize,
        p=c(0.05,0.5,0.95),q=quantile(cases,prob=p,names=FALSE)) %>%
  mutate(p=mapvalues(p,from=c(0.05,0.5,0.95),to=c("lo","med","hi")),
         data=mapvalues(data,from=c(TRUE,FALSE),to=c("data","simulation"))) %>%
  dcast(time+data~p,value.var='q') %>%
  ggplot(aes(x=time,y=med,color=data,fill=data,ymin=lo,ymax=hi))+
  geom_ribbon(alpha=0.2)

The parameters are constrained to be positive, and some of them are constrained to lie between \(0\) and \(1\). We can turn the likelihood maximization problem into an unconstrained maximization problem by transforming the parameters. The following Csnippets implement such a transformation and its inverse.

toEst <- Csnippet("
  Tmu = log(mu);
  Tsigma = log(sigma);
  Tgamma = log(gamma);
  Talpha = log(alpha);
  Tiota = log(iota);
  Trho = logit(rho);
  Tcohort = logit(cohort);
  Tamplitude = logit(amplitude);
  TsigmaSE = log(sigmaSE);
  Tpsi = log(psi);
  TR0 = log(R0);
  to_log_barycentric (&TS_0, &S_0, 4);
")

fromEst <- Csnippet("
  Tmu = exp(mu);
  Tsigma = exp(sigma);
  Tgamma = exp(gamma);
  Talpha = exp(alpha);
  Tiota = exp(iota);
  Trho = expit(rho);
  Tcohort = expit(cohort);
  Tamplitude = expit(amplitude);
  TsigmaSE = exp(sigmaSE);
  Tpsi = exp(psi);
  TR0 = exp(R0);
  from_log_barycentric (&TS_0, &S_0, 4);
")


m1 <- pomp(m1,toEstimationScale=toEst,
           fromEstimationScale=fromEst,
           statenames=c("S","E","I","R","C","W"),
           paramnames=c("R0","mu","sigma","gamma","alpha","iota",
                        "rho","sigmaSE","psi","cohort","amplitude",
                        "S_0","E_0","I_0","R_0"))

References

He, D., E. L. Ionides, and A. A. King. 2010. Plug-and-play inference for disease dynamics: Measles in large and small populations as a case study. Journal of The Royal Society Interface 7:271–283.